% In this practical we will work with a single subject's EEG dat and perform an ERF
% analysis in sensor space. 
%

%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

clear all; close all;

global OSLDIR;


tilde='/Users/dantemant/Documents/meeg';
osldir=[tilde '/osl1.2.beta.16dante'];    

addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS

mri_dir='/Users/dantemant/Documents/meeg/meeg_data/mr_images';

testdir=[tilde '/meeg_data/meeg_signals/case_1084/111202'];

workingdir=[testdir '/preprocessing']; % this is the directory the SPM files will be stored in

cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in

clear spm_files_continuous spm_files_epoched;

% set up a list of SPM MEEG object file names (we only have one here)
spm_files_continuous{1}=[workingdir '/CSMCAdfMspm8_eeg1.mat'];
spm_files_epoched{1}=[workingdir '/eCSMCAdfMspm8_eeg1.mat'];


% defining experimental conditions and contrasts to be calculated

% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to croat_settingse the (num_regressors x num_trials) design matrix:


cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in

% Set up the list of subjects and their structural scans for the analysis 
% Currently there is only 1 subject.

% set up a list of mris
structural_files{1}=[mri_dir '/struct_S01.nii'];

% set up a list of SPM MEEG object file names (we only have one here)

%%%%%%%%%%%%%%%%%%%
%% DO REGISTRATION AND RUN FORWARD MODEL BASED ON STRUCTURAL SCANS
% Before running the beamformer we need to compute the forward model for
% each subject based on their structural scan.
% Make sure you check the results look reasonable!

for i=1:length(spm_files_continuous),

    %D=spm_eeg_load(spm_files_epoched{i});
    %fidnew=D.fiducials;
    %fidnew.fid.label

    S2=[];

    S2.fid_label.nasion='Nasion';        
    S2.fid_label.lpa='LPA';
    S2.fid_label.rpa='RPA';

    S2.D = spm_files_continuous{i};    % requires .mat extension    
    %S2.D = spm_files_epoched{i};    % requires .mat extension    
    
    S2.mri=structural_files{i}; % set S2.mri=''; if there is no structural available        
    S2.useheadshape=1;

    %S2.forward_meg='Single Shell';
    %S2.forward_meg='MEG Local Spheres';

    D=osl_forward_model(S2);

    %D=osl_neuromag_grad_baseline_correction(S2.D,'vector_view');

end;

%spm_eeg_inv_checkmeshes(D);
spm_eeg_inv_checkdatareg(D);
%spm_eeg_inv_checkforward(D, 1);



%% defining experimental conditions and contrasts to be calculated

% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to croat_settingse the (num_regressors x num_trials) design matrix:

conditions=[{'9'},{'25'},{'41'},{'1'}, {'17'},{'33'}];

% conditions=[{0},{1},{2},{3},{4},{5}, ...
%     {6},{7},{8},{9},{10},{11},{12}, ...
%     {13},{14},{15},...
%     {20},{21},{22},{23},{24},{25}, ...
%     {26},{27},{28},{29},{30},{31},{32}, ...
%     {33},{34},{35},...
%     {40},{41},{42},{43},{44},{45}, ...
%     {46},{47},{48},{49},{50},{51},{52}, ...
%     {53},{54},{55}];



Xsummary={};
cond_list=[1:6];
for i=1:length(cond_list),
    Xsummary{i}=zeros(size(conditions,2),1);
    Xsummary{i}(cond_list(i))=1;
end

contrast={};

contrast{1}=zeros(length(cond_list),1);
contrast{1}([1 2 3])=1; % all targets
contrast{2}=zeros(length(cond_list),1);
contrast{2}([4 5 6])=1; % away non-targets


% setup oat
oat=[];
%oat.source_recon.D_continuous{1}=spm_files_continuous{1};
oat.source_recon.conditions=conditions;
oat.source_recon.D_epoched{1}=spm_files_epoched{1}; % this is passed in so that the bad trials and bad channels can be read out
oat.source_recon.freq_range=[]; % frequency range in Hz
oat.source_recon.time_range=[-0.1 0.5];
oat.source_recon.modalities={'EEG'};
oat.source_recon.method='beamform';
oat.source_recon.gridstep=12; % in mm, using a lower resolution here than you would normally, for computational speed
oat.source_recon.mri=structural_files;
oat.source_recon.dirname=[spm_files_continuous{1} '_wideband'];
%oat.source_recon.forward_meg='MEG Local Spheres';
oat.source_recon.forward_meg='Single Shell';
oat.source_recon.work_in_pca_subspace=0;
oat.source_recon.pca_dim=50;
%oat.source_recon.hmm_num_states=-1;
%oat.source_recon.hmm_num_starts=1;
oat.first_level.design_matrix_summary=Xsummary;

for i=1:length(contrast)
    oat.first_level.contrast{i}=contrast{i};
    oat.first_level.contrast_name{i}=['C' num2str(i)];
end
oat.first_level.cope_type='acope';

oat = osl_check_oat(oat);  % function changed by DM


%% RUN THE OAT:

oat.to_do=[1 1 0 0];
oat.first_level.name=['wholebrain_first_level'];
oat.do_plots=1;
oat = osl_run_oat(oat);

% load GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% view the GLM design matrix
figure;imagesc(stats.x);title('GLM design matrix');xlabel('regressor no.');ylabel('trial no.');

% results = osl_get_recon_timecourses( osl_load_oat_results(oat,oat.source_recon.results_fnames{1});, mask_fname )

%% OUTPUT SUBJECT'S NIFTII FILES
% Having run the GLM on our source space data, we would like to inspect the
% results for our single subject. 
% We can do this by saving the contrast of parameter estimates (COPEs) and 
% t-statistics for each of our contrasts to NIFTI images.

S2=[];
S2.time.reduce_time=1;
S2.resample_method='fft_ds';
S2.time_ds=2;
S2.time_range=[-0.1 0.3];
S2.data=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
[stats] = osl_reduce_data_to_visualize(S2);

S2=[];
S2.oat=oat;
S2.stats=stats;
S2.resamp_gridstep=12;
S2.first_level_contrasts=[1]; % list of first level contrasts to output
[statsdir,times,count]=osl_save_nii_stats(S2);

%% VIEW NIFTII RESULTS IN FSLVIEW
% We can now view the nifti images containing our GLM results in FSL, here
% we are runnin fslview from the matlab command line, but you do not need
% to - you can run it from the UNIX command line instead.


con=1;
runcmd([getenv('FSLDIR') '/bin/fslview ' [statsdir '/tstat' num2str(con) '_12mm']  ' &']);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
